Abstract

This document is for the analysis of data collected in Southeast Alaska to assess variation in food web structure in streams across a gradient of hydrologic and physio-chemical conditions. This is a workind document and will be backed up on github at the following location: https://github.com/mdunkle/SEAK-Analysis.git

Fish Data

First thing I want to do is figure out the total biomass of fish captured in each sampling event.

## # A tibble: 455 x 19
## # Groups:   ReachID, Date [15]
##     ReachID      StreamID       Date Species Length Weight  Temp Wwidth
##      <fctr>        <fctr>     <date>  <fctr>  <int>  <dbl> <int>  <dbl>
##  1 MT_Upper Montana Creek 2017-07-21    Coho     39   0.60    12    8.9
##  2 MT_Upper Montana Creek 2017-07-21    Coho     45   1.30    12    8.9
##  3 MT_Upper Montana Creek 2017-07-21    Coho     44   1.10    12    8.9
##  4 MT_Upper Montana Creek 2017-07-21    Coho     42   1.00    12    8.9
##  5 MT_Upper Montana Creek 2017-07-21    Coho     37   0.50    12    8.9
##  6 MT_Upper Montana Creek 2017-07-21    Coho     79   0.56    12    8.9
##  7 MT_Upper Montana Creek 2017-07-21    Coho     84   0.73    12    8.9
##  8 MT_Upper Montana Creek 2017-07-21    Coho     44   1.00    12    8.9
##  9 MT_Upper Montana Creek 2017-07-21    Coho     38   0.60    12    8.9
## 10 MT_Upper Montana Creek 2017-07-21    Coho     54   1.90    12    8.9
## # ... with 445 more rows, and 11 more variables: BFWidth <dbl>,
## #   Depth <dbl>, StreamType <fctr>, SamplingEvent <fctr>, X <lgl>,
## #   Species.1 <lgl>, CPUE <lgl>, StreamID.1 <lgl>, ReachID.1 <lgl>,
## #   number_fish <int>, biomass <dbl>

Plotting the length/weight relationship for Coho, Cutthroat, and DVs at the MT/McginConf

Here I’m visualizing lengths by species and reach

Coho Data

Next, I’m visualizing just the Dolly Varden data and showing the distribution as boxplots.

Dolly Varden Data

Next, I’m visualizing just the Dolly Varden data and showing the distribution as a boxplots.

Length Distribution by Species and Sampling Event

Here, I’m using density ridge plots to show the distribution of lengths of each species by the two sampling events. The second plot shows the same data, but bracketed by stream type (glacial, clearwater (cw), brownwater (bw), and a combination of two (e.g. (CW/BW)) or all three (combined)). These distinctions are somewhat arbitrary, as most systems in SEAK show a combination of all three stream types and no one system is purely one type. It may be more useful to represent each of these as bins of conditions along a gradient of turbidity, flow variability, or DOC concentration.

cohodv_data = SEAK %>% filter(Species == c("DollyVarden","Coho")) %>% group_by(StreamID, SamplingEvent, Species)

cohodv=ggplot(SEAK, aes(x=Length, y=StreamID, fill=Species))+
  geom_point(aes(color=Species), alpha=.5)+
  geom_density_ridges(alpha=.5,scale=1,size=.00001, rel_min_height=0)+facet_wrap(~SamplingEvent)+
  theme_cowplot()+theme(strip.background = element_blank(), legend.title = element_blank())
cohodv

SEAK_MeanWeight_Summary = SEAK %>% group_by(ReachID, StreamID, SamplingEvent, Species) %>% dplyr::summarise(mean=mean(Weight), na.rm=T)

ggplot(SEAK_MeanWeight_Summary, aes(x=SamplingEvent, y = mean, color=Species, group=Species))+geom_point()+xlab(NULL)+ylab("Mean Weight")+
  facet_grid(~StreamID, switch="both")+
  theme(axis.text.x=element_blank(), axis.ticks.x = element_blank(),
        strip.background = element_blank(), strip.placement = "outside", strip.text.x = element_text(angle=90, vjust=1),
        panel.spacing.x=unit(0, "lines"))+
  geom_path(aes(x=SamplingEvent, y = mean, group=Species), arrow = arrow(type="closed"), size=1, position = position_dodge(width=1), alpha=.5)
## Warning: Removed 16 rows containing missing values (geom_point).
## Warning: Removed 10 rows containing missing values (geom_path).
## geom_path: Each group consists of only one observation. Do you need to
## adjust the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to
## adjust the group aesthetic?

Catch Per Unit Effort Across the Streams

Here, I load the data for catch per unit effort stored in the file SEAK_CPUE.csv. The associated figure shows the average catch across sites within a stream system during July and September sampling with an arrow showing the trend in catch rate between the two events. Important to note, Peterson Creek had the highest catch rate (and likely density) of any site in July and almost no juvnile salmonids in September. This could be due to hypoxic conditions or disturbance caused by adult salmon spawning in high densities in this system.

SEAK_CPUE=read.csv("SEAK_CPUE.csv", head=T)
SEAK_CPUE_Summary = SEAK_CPUE %>% group_by(StreamID, SamplingEvent, Species) %>% dplyr::summarise(mean=mean(CPUE))

print(head(SEAK_CPUE_Summary))
## # A tibble: 6 x 4
## # Groups:   StreamID, SamplingEvent [2]
##     StreamID SamplingEvent           Species  mean
##       <fctr>        <fctr>            <fctr> <dbl>
## 1 Fish Creek          July           Chinook     1
## 2 Fish Creek          July CoastRangeSculpin     3
## 3 Fish Creek          July              Coho     9
## 4 Fish Creek          July         Cutthroat     0
## 5 Fish Creek          July       DollyVarden     2
## 6 Fish Creek     September           Chinook     0
catch_trends=ggplot(SEAK_CPUE_Summary, aes(x=SamplingEvent, y = mean, color=Species, group=SamplingEvent))+
  #geom_point()+
  xlab(NULL)+ylab("Mean CPUE")+
  facet_grid(~StreamID, switch="both")+
  theme(axis.text.x=element_blank(), axis.ticks.x = element_blank(),legend.title=element_blank(),
        strip.background = element_blank(), strip.placement = "outside", strip.text.x = element_text(angle=90, vjust=1),
        panel.spacing.x=unit(0, "lines"))+scale_y_log10()+
  geom_path(aes(x=SamplingEvent, y = mean, group=Species), arrow = arrow(type="closed"), size=1, position = position_dodge(width=1), alpha=.5)
catch_trends

save_plot("2017_CatchTrends.pdf", catch_trends, base_width = 8, base_height = 6)

Periphyton Data

In this section, I visualize and analyze the data from periphyton samples taken during the 2017 pilot study.

Auxilliary Data

UAS Data from Jason Fellman showing the variables measured in 2006-2007.

## `geom_smooth()` using method = 'loess'
## Warning: Removed 7 rows containing non-finite values (stat_smooth).
## `geom_smooth()` using method = 'loess'
## Warning: Removed 6 rows containing non-finite values (stat_smooth).
## `geom_smooth()` using method = 'loess'
## Warning: Removed 15 rows containing non-finite values (stat_smooth).

Modeling Simlulations

This was a practice for extracting data from Stella simulations and plotting in ggplot.

Making Maps for the proposal

Here I’m doing some spatial work to make sampling location maps. For some reason, this code is no longer working and needs some attention.

r, echo=F, message=F, warning = F} #library(“tmap”) library(“tmaptools”) library(“sf”) library(“leaflet”) library(“rio”) library(“maptools”)

wbdhu14 = “WBDHU14.shp” wbdhu14geo=read_shape(file=wbdhu14, as.sf=T) qtm(wbdhu14geo)

wbdhu8 = “WBDHU8.shp” wbdhu8geo = read_shape(file=wbdhu8, as.sf=T) qtm(wbdhu8geo)

wbdhu12 = “WBDHU12.shp” wbdhu12geo = read_shape(file=wbdhu12, as.sf=T) qtm(wbdhu12geo, “NAME” )

wbdhu10=“WBDHU10.shp” wbdhu10geo=read_shape(file=wbdhu10, as.sf = T) qtm(wbdhu10geo)

nhdflowline=“NHDFlowline.shp” nhdflowlinegeo=read_shape(file=nhdflowline, as.sf=T) qtm(nhdflowlinegeo)

alaskacoastline=“alaska_63360_py.shp” alaskacoastlinegeo=read_shape(file=alaskacoastline) pcreek=“Peterson Creek” wbdhu14geo=cbind(pcreek, wbdhu14geo) seakpoints=“SEAK-Sampling-Locations.shp” seakpointsgeo=read_shape(file=seakpoints) juneauwater=“tl_2015_02110_areawater.shp” juneauwatergeo=read_shape(file=juneauwater) head(juneauwatergeo, 50) mendenhallglacier=subset(juneauwatergeo, FULLNAME==“Mendenhall Glacier”) juneauicefield=subset(juneauwatergeo, FULLNAME==“Juneau Icefield”) akglaciers = “01_rgi60_Alaska.shp” akglaciersgeo=read_shape(file=akglaciers)

tm=tm_shape(wbdhu12geo, ylim=c(58.35, 58.55), xlim=c(-134.8, -134.38))+ tm_fill(“white”)+tm_shape(wbdhu12geo[wbdhu12geo$NAME==“Lynn Canal”,])+tm_fill(“gray20”)+ #tm_shape(wbdhu12geo[wbdhu12geo$NAME==“Eagle River”,])+tm_fill(“white”, lwd=2)+tm_text(text=“NAME”)+tm_borders(“black”)+ # tm_shape(wbdhu12geo[wbdhu12geo$NAME==“Herbert Glacier”,])+tm_fill(“white”, lwd=2)+tm_borders(“black”)+ #tm_shape(wbdhu12geo[wbdhu12geo$NAME==“Mendenhall Glacier”,])+tm_fill(“white”, lwd=2)+tm_borders(“black”)+ tm_shape(wbdhu14geo[wbdhu14geo$NAME==“Fritz Cove-Frontal Lynn Canal”,])+tm_fill(“gray20”, lwd=2)+tm_borders(“gray20”)+ tm_shape(wbdhu14geo[wbdhu14geo$NAME==“Tee Creek-Frontal Lynn Canal”,])+tm_fill(“gray20”, lwd=2)+tm_borders(“gray20”)+ tm_shape(wbdhu14geo[wbdhu14geo$NAME==“19010301060612-Mendenhall Peninsula”,])+tm_fill(“gray20”, lwd=2)+ tm_shape(alaskacoastlinegeo)+tm_fill(“white”)+ tm_shape(juneauicefield)+tm_fill(“dodgerblue”)+#tm_text(“FULLNAME”)+ tm_shape(mendenhallglacier)+tm_fill(“dodgerblue”)+#tm_text(“FULLNAME”)+ tm_shape(akglaciersgeo)+tm_fill(col=“dodgerblue”)+ tm_shape(wbdhu14geo[wbdhu14geo$NAME==“Mendenhall River”,])+tm_fill(“gray”, lwd=2, alpha=.5)+tm_text(text=“NAME”, xmod=6, ymod=4)+tm_borders(“black”)+ tm_shape(wbdhu14geo[wbdhu14geo$NAME==“Steep Creek”,])+tm_fill(“gray”, lwd=2, alpha=.5)+tm_text(text=“NAME”, xmod=3, ymod=2)+tm_borders(“black”)+ tm_shape(wbdhu12geo[wbdhu12geo$NAME==“Herbert River”,])+tm_fill(“gray”, lwd=2, alpha=.5)+tm_text(text=“NAME”)+tm_borders(“black”)+ tm_shape(wbdhu14geo[wbdhu14geo$NAME==“Peterson Lake-Peterson Creek”,])+tm_fill(“gray”, lwd=2, alpha=.5)+tm_text(text=“pcreek”, xmod=1, ymod=-1)+tm_borders(“black”)+ tm_shape(wbdhu12geo[wbdhu12geo$NAME==“Montana Creek”,])+tm_fill(“gray”, lwd=2, alpha=.5)+tm_text(text=“NAME”, xmod=0,ymod=-6)+tm_borders(“black”)+ #tm_shape(wbdhu14geo[wbdhu14geo$NAME==“Windfall Creek”,])+tm_fill(“white”,lwd=2)+tm_text(text=“NAME”)+tm_borders(“black”)+ tm_shape(wbdhu14geo[wbdhu14geo$NAME==“McGinnis Creek”,])+tm_fill(“gray”,lwd=2, alpha=.0)+tm_text(text=“NAME”, xmod=7.5, ymod=2)+tm_borders(“black”)+ #tm_grid(projection=“longlat”, labels.size=.5, alpha=.75)+ tm_shape(nhdflowlinegeo)+tm_lines(“black”, alpha=.5)+ tm_shape(seakpointsgeo)+tm_symbols(col=“black”)+ tm_scale_bar()+tm_style_classic() tm

save_tmap(tm, “AKFieldsites.png”, width=3000, height=3000)

# function to obtain US county shape get_US_county_2010_shape <- function() { dir <- tempdir() download.file(“http://www2.census.gov/geo/tiger/GENZ2010/gz_2010_us_050_00_20m.zip”, destfile = file.path(dir, “gz_2010_us_050_00_20m.zip”)) unzip(file.path(dir, “gz_2010_us_050_00_20m.zip”), exdir = dir) read_shape(file.path(dir, “gz_2010_us_050_00_20m.shp”)) }

obtain US county shape

US <- get_US_county_2010_shape()

split shape

US_cont <- US[!(US$STATE %in% c(“02”,“15”,“72”)),]
US_AK <- US[US$STATE == “02”, ] US_HI <- US[US$STATE == “15”,]

create state boundaries

US_states <- unionSpatialPolygons(US_cont, IDs=US_cont$STATE) tmap_mode(“plot”)

m_AK <- tm_shape(US_AK, projection = 3338) + tm_polygons(border.col = “grey50”, border.alpha = .5, breaks = seq(10, 50, by = 5)) + tm_layout(“Alaska”, legend.show = FALSE, bg.color = NA, title.size = 0.8, frame = FALSE) print(m_AK) ```